autocorrelation
Model 3: SEM CM without respecifications
Results
## Soil_property ME RMSE SS
## 1 CEC.A 0.00775128089893602 5.68793208043344 5047.00113085333
## 2 CEC.B -0.00398182488814064 6.26155391338747 6116.30095600024
## 3 CEC.C 0.00229007311875548 6.77752944249498 7165.84523364626
## 4 OC.A 0.000480842199394078 0.629623927110415 61.8425011760308
## 5 OC.B 0.000916586572794572 0.316970003100663 15.673317327039
## 6 OC.C 0.000757381666383179 0.180989934293459 5.11014758522591
## 7 clay.A 0.00387041702804282 8.34581735151287 10865.8160933109
## 8 clay.B -0.0125981787909463 9.39142974125339 13759.0366032442
## 9 clay.C -0.000429563387851167 10.0827907785092 15859.3765017776
## mean_theta median_th R2
## 1 0.9450710 0.1988536 0.05396970
## 2 0.7638961 0.2693472 0.15050286
## 3 0.7684111 0.3220226 0.16161086
## 4 1.0386352 0.3081390 0.06873298
## 5 0.9604539 0.2237897 0.02997153
## 6 1.4402485 0.1209572 0.04206266
## 7 0.8671333 0.1784460 0.11736397
## 8 0.8878465 0.3450671 0.21954429
## 9 1.1430615 0.2243374 0.15005788## Soil_property ME RMSE r2
## 1 CECo 0.0020198430 6.258182 0.2025248
## 3 OCo 0.0007182701 0.420180 0.6836467
## 5 clayo -0.0030524417 9.300794 0.2659860Spatial autocorrelation analysis
### Autocorrelation in SEM residuals ####
backup <- Res
Res[2:10] <- Res[,2:10]-Res[,11:19]
Res <- Res[,1:10]
names(Res)[1] <- "id.p"
R <- merge(Res, unique(d[,c(1,20,21)]), by.x = "id.p", by.y = "idp", all.x=T)
library(sp)
library(gstat)
# R as geo
coordinates(R) <- ~X+Y
#define crs
wgs84 <- CRS("+init=epsg:4326")
#UTM14N <- CRS("+init=epsg:32614")
modis <- CRS("+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs")
NAD83.KS.N <- CRS("+init=epsg:2796")
proj4string(R) <- NAD83.KS.N
R <- spTransform(R, wgs84)
mapview::viewExtent(R)
library(mapview)## Loading required package: leafletmapview(R, burst = TRUE)###############################
# Variograms
# independent of device size
g <- list()
vg <- list()
vgm <- list()
# CEC
for(i in 2:4){
g[[i-1]] <- gstat(id = c(names(R@data)[i]),
formula = formula(paste0(names(R@data)[i],"~1")),
data = R)
vg[[i-1]] <- variogram(g[[i-1]])
# vg = variogram(g, width = 20000, cutoff = 600000)
# vg = variogram(g, boundaries = c(1E4,2E4,4E4,7E4,1E5,2E5,4E5,7E5,1E6))
# print(plot(vg, plot.numbers = TRUE, main = names(R@data)[i]))
# # choose initial variogram model and plot:
vgm[[i-1]] <- vgm(nugget = var(R@data[,i]),
psill=var(R@data[,i])*2,
range=100000,
model = "Sph")
#plot(vg, vgm, main = names(R@data)[i])
#
# # fit variogram model:
vgm[[i-1]] <- fit.variogram(vg[[i-1]], vgm[[i-1]])
print(names(R@data)[i])
print(vgm[[i-1]])
print(attr(vgm[[i-1]], "SSErr"))
}## Warning in fit.variogram(vg[[i - 1]], vgm[[i - 1]]): No convergence after
## 200 iterations: try different initial values?## [1] "CEC.A"
## model psill range
## 1 Nug 32.56124 0.0
## 2 Sph 65.12263 100000.2
## [1] 487.1739## Warning in fit.variogram(vg[[i - 1]], vgm[[i - 1]]): No convergence after
## 200 iterations: try different initial values?## [1] "CEC.B"
## model psill range
## 1 Nug 39.45999 0.0
## 2 Sph 78.91982 99999.8
## [1] 1681.899## Warning in fit.variogram(vg[[i - 1]], vgm[[i - 1]]): No convergence after
## 200 iterations: try different initial values?## [1] "CEC.C"
## model psill range
## 1 Nug 46.23125 0.00
## 2 Sph 92.46235 99999.83
## [1] 4387.061# OC
for(i in 5:7){
g[[i-1]] <- gstat(id = c(names(R@data)[i]),
formula = formula(paste0(names(R@data)[i],"~1")),
data = R)
vg[[i-1]] <- variogram(g[[i-1]])
# vg = variogram(g, width = 20000, cutoff = 600000)
# vg = variogram(g, boundaries = c(1E4,2E4,4E4,7E4,1E5,2E5,4E5,7E5,1E6))
# print(plot(vg, plot.numbers = TRUE, main = names(R@data)[i]))
# # choose initial variogram model and plot:
vgm[[i-1]] <- vgm(nugget = var(R@data[,i]),
psill= var(R@data[,i])*2,
range=100000,
model = "Sph")
#plot(vg, vgm, main = names(R@data)[i])
#
# # fit variogram model:
vgm[[i-1]] <- fit.variogram(vg[[i-1]], vgm[[i-1]])#, fit.method = 4)
print(names(R@data)[i])
print(vgm[[i-1]])
print(attr(vgm[[i-1]], "SSErr"))
}## Warning in fit.variogram(vg[[i - 1]], vgm[[i - 1]]): No convergence after
## 200 iterations: try different initial values?## [1] "OC.A"
## model psill range
## 1 Nug 0.3989836 0.0
## 2 Sph 0.7979665 99999.9
## [1] 0.1327332## Warning in fit.variogram(vg[[i - 1]], vgm[[i - 1]]): No convergence after
## 200 iterations: try different initial values?## [1] "OC.B"
## model psill range
## 1 Nug 0.1011173 0.0
## 2 Sph 0.2022348 100000.1
## [1] 0.008639608## Warning in fit.variogram(vg[[i - 1]], vgm[[i - 1]]): No convergence after
## 200 iterations: try different initial values?## [1] "OC.C"
## model psill range
## 1 Nug 0.03296812 0.0
## 2 Sph 0.06593627 100000.1
## [1] 0.0009133788# Clay
for(i in 8:10){
g[[i-1]] <- gstat(id = c(names(R@data)[i]),
formula = formula(paste0(names(R@data)[i],"~1")),
data = R)
vg[[i-1]] <- variogram(g[[i-1]])
# vg = variogram(g, width = 20000, cutoff = 600000)
# vg = variogram(g, boundaries = c(1E4,2E4,4E4,7E4,1E5,2E5,4E5,7E5,1E6))
# print(plot(vg, plot.numbers = TRUE, main = names(R@data)[i]))
# # choose initial variogram model and plot:
vgm[[i-1]] <- vgm(nugget = var(R@data[,i]),
psill=var(R@data[,i])*2,
range=5E4,
model = "Sph")
#plot(vg, vgm, main = names(R@data)[i])
#
# # fit variogram model:
vgm[[i-1]] <- fit.variogram(vg[[i-1]], vgm[[i-1]], fit.method = 2)
print(names(R@data)[i])
print(vgm[[i-1]])
print(attr(vgm[[i-1]], "SSErr"))
}## Warning in fit.variogram(vg[[i - 1]], vgm[[i - 1]], fit.method = 2): No
## convergence after 200 iterations: try different initial values?## [1] "clay.A"
## model psill range
## 1 Nug 70.10202 0.0
## 2 Sph 140.20911 50001.8
## [1] 336.8401## Warning in fit.variogram(vg[[i - 1]], vgm[[i - 1]], fit.method = 2): No
## convergence after 200 iterations: try different initial values?## [1] "clay.B"
## model psill range
## 1 Nug 88.76782 0.00
## 2 Sph 177.53011 49998.44
## [1] 299.7812## Warning in fit.variogram(vg[[i - 1]], vgm[[i - 1]], fit.method = 2): No
## convergence after 200 iterations: try different initial values?## [1] "clay.C"
## model psill range
## 1 Nug 102.3186 0.00
## 2 Sph 204.6411 50000.97
## [1] 621.602# Three graphs (soil properties)
par(mfrow = c(1, 3), pty = "s")
### CEC
## A
plot(variogramLine(vgm[[1]], maxdist=2E5),
type="l", lwd=2,col="#AA0000", main="CEC",
xlab = "Distance", ylab = "Semivariance", cex.lab = 1.3, ylim=c(0,50))
points(gamma ~ dist, vg[[1]], col="#770000")
#legend(x= "topleft",legend = "SSErr", bty = "n")
#legend(x= 12,legend = round(attr(vgm[[1]], "SSErr"),3), text.col ="#AA0000",
# bty = "n")
## B
lines(variogramLine(vgm[[2]], maxdist=2E5), lwd=2, col="#00AA00")
points(gamma ~ dist, vg[[2]], col="#007700")
## C
lines(variogramLine(vgm[[3]], maxdist=2E5), lwd=2, col="#0000AA")
points(gamma ~ dist, vg[[3]], col="#000077")
### OC
## A
plot(variogramLine(vgm[[4]], maxdist=2E5), type="l", lwd=2,col="#AA0000",
main="OC",
xlab = "Distance", ylab = "Semivariance", cex.lab = 1.3, ylim=c(0,0.5))
points(gamma ~ dist, vg[[4]], col="#770000")
## B
lines(variogramLine(vgm[[5]], maxdist=2E5), lwd=2, col="#00AA00")
points(gamma ~ dist, vg[[5]], col="#007700")
## C
lines(variogramLine(vgm[[6]], maxdist=2E5), lwd=2, col="#0000AA")
points(gamma ~ dist, vg[[6]], col="#000077")
### Clay
## A
plot(variogramLine(vgm[[7]], maxdist=2E5), type="l", lwd=2,col="#AA0000",
main="Clay",
xlab = "Distance", ylab = "Semivariance", cex.lab = 1.3, ylim=c(0,140))
points(gamma ~ dist, vg[[7]], col="#770000")
## B
lines(variogramLine(vgm[[8]], maxdist=2E5), lwd=2, col="#00AA00")
points(gamma ~ dist, vg[[8]], col="#007700")
## C
lines(variogramLine(vgm[[9]], maxdist=2E5), lwd=2, col="#0000AA")
points(gamma ~ dist, vg[[9]], col="#000077")
#############################################################################
par(mfrow = c(1, 1))
title("Semivariograms of residuals",
sub="red: A horizon, green: B horizon, blue: C horizon")
# vgm[[4]] <- vgm(nugget = 0.1,
# psill=0.1,
# range=5E4,
# model = "Wav")
# plot(vg[[4]], vgm[[4]],ylim=c(0,0.5))
# vgm[[4]] <- fit.variogram(vg[[4]], vgm[[4]], fit.method = 2)
# plot(vg[[4]], vgm[[4]],ylim=c(0,0.5))
# attr(vgm[[4]], "SSErr")
# vgm[[4]]